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Molecular dynamics computer simulations are used to study the aging dynamics of Si02 (mod- 
eled by the BKS model). Starting from fully equilibrated configurations at high temperatures 
Ti e {5000 K, 3760 K}, the system is quenched to lower temperatures Tf € {2500 K, 2750 K, 3000 K, 
3250 K} and observed after a waiting time tw Since the simulation runs are long enough to reach 
equilibrium at Tf , we are able to study the transition from out-of-equilibrium to equilibrium dynam- 
ics. We present results for the partial structure factors, for the generalized incoherent intermediate 
scattering function Cq{tv,,t„+t), and for the mean square displacement Ar^(tw, iw+i). We conclude 
that there are three different regions: (I) At very short waiting times, Cq{t^,t„ + 1) decays very 
fast without forming a plateau. Similarly Ar'^(t„,iw +t) increases without forming a plateau. (II) 
With increasing tw a plateau develops in Cq{ty„,t„-\-t) and Ar^(tw, tw + t). For intermediate waiting 
times the plateau height is independent of tw and Ti. Time superposition applies, i.e. Cq — Cq{t/t^'^) 
where t^"^ = t^'^{t^) is a waiting time dependent decay time. Furthermore Cq = C{q,t^,t^ + t) 
scales as Cq = C{q, z{t^,t) where 2 is a function of tw and t only, i.e. independent of q. (Ill) At 
large tw the system reaches equilibrium, i.e. Cq{t„,t^ + t) and Ar^(fw,tw + t) are independent of 
tw and Ti. For Cq{t^,t^ -\-t) we find that the time superposition of intermediate waiting times (II) 
includes the equilibrium curve (III). 

PACS numbers: 61.20.Lc, 61.20.Ja, 64.70.ph, 02.70.Ns, 61.43.Fs 
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I. INTRODUCTION 

When a glass-forming liquid is quenched from an 
equilibrium state at a high temperature Ti to a non- 
equilibrium state at a lower temperature Tf, "aging pro- 
cesses" set in. Provided that crystallization plays no role 
at Tf (e.g. due to very low crystal nucleation rates), the 
transition to a final (metastable) equilibrium state occurs 
on a time scale that corresponds to the typical equilib- 
rium relaxation time Teq of the (supercooled) liquid at 
Tf. The dynamics of the system depends on the wait- 
ing time which is the time elapsed after the tempera- 
ture quench. If Tcq exceeds the waiting time then the 
system is observed in a transient non-equilibrium state 
which corresponds to a glass for <^ Tcq. During the 
aging process, i.e. for < Toq, thermodynamic proper- 
ties such as volume and energy are changing and time 
translation invariance does not hold: correlation func- 
tions at time ty^ +t and the time origin at do depend 
not only on the time difference t but also on the waiting 
time tw 

Recently this aging process has been investigated ex- 
tensively with experiments [IHll, theoretically [3-[l| and 
with computer simulations. For a more complete sum- 
mary of previous results we refer the reader to the ref- 
erences [13, [HI and references therein. Computer sim- 
ulation studies most similar to the work presented here 
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are on attractive colloidal systems |12h15| . on the Kob- 
Andersen Lennard- Jones (KALJ) mixture [T6l - [25j . and 
silica (8102) [25l - l28j . In the case of silica the interpreta- 
tion of the results is less clear than for the KALJ mixture; 
e.g. different findings [1^ H^l have been reported on the 
violation of the fluctuation-dissipation regime during ag- 
ing Thus, it remains open whether silica, as the 
prototype of a glass-forming system forming a tetrahe- 
dral network structure, exhibits a different aging dynam- 
ics than, e.g. the KALJ model, where the structure is 
similar to that of a closed-packed hard-sphere structure. 

Recent simulation studies on amorphous silica (25l - [39| 
have widely used the BKS potential [i^ to model the 
interactions between the atoms. Although it is a simple 
pair potential, it reproduces various static and dynamic 
properties of amorphous silica very well. For BKS sil- 
ica, the self-diffusion constants Da {a = Si, O) show two 
different temperature regimes: At high temperatures. 
Da decays according to a power law, as predicted by 
the mode coupling theory (MCT) of the glass transition 
(note, however, that also other interpretations have been 
assigned to this high temperature regime). At low tem- 
perature. Da as well as the shear viscosity 77 exhibit an 
Arrhenius behavior with an activation ener gy o f the order 
of 5 eV, in agreement with experiment (see [3C| and refer- 
ences therein). The temperature at which the crossover 
between both regimes occurs is at Tc « 3300 K, corre- 
sponding to the critical MCT temperature of BKS silica. 
Previous studies of the aging dynamics of BKS silica [25|- 
[27! were performed in two steps. First, the system was 
fully equilibrated at a temperature Ti > T^. Then, the 
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system was quenched to a low temperature Tf < Tc, fol- 
lowed by the production runs. Wahlen and Rieger 
analyze time-dependent correlation functions at different 
waiting times and Berthier [25'] and Scala et al. ^2l\ 
study the generalized fluctuation dissipation relation and 
the energy landscape [l^l- All three studies p5l - l27l | in- 
vestigate the early stages of the aging dynamics, i.e. the 
dynamics was explored on time scales that were much 
smaller than the equilibrium relaxation time Tcq at the 
temperature T(. 

In this work, we also consider quenches in BKS sil- 
ica from a high temperature Tj to a low temperature Tf . 
Different from previous simulation studies, we aim at elu- 
cidating the full transient dynamics at Tf from the initial 
state at = to equilibrium. To this end, temperatures 
Tf are chosen such that the system can be fully equili- 
brated on the typical time span of the MD simulation. 
Note that the considered temperatures Tf G {2500 K, 
2750 K, 3000 K, 3250 K} are below the critical MCT tem- 
peratures Tc. Thus, we have access to the full aging dy- 
namics in the experimentally relevant Arrhenius temper- 
ature regime that we have mentioned above. 

The analysis of time-dependent density correlation 
functions Cq {t„ ,t„ + t) (with q the wavenumber) and the 
mean square displacement Ar^ (iw ,tw + t) reveal three 
different regimes of waiting times In the case of 
Cq{t„, t„ + t) (and similarly for Ar^(iw, + *)) at early 
t„, a rapid decay to zero is seen, without forming a 
plateau at intermediate times. Then, for larger values 
of a plateau is formed. The height of this plateau 
grows with waiting time and becomes more pronounced, 
before in the final regime the plateau height is indepen- 
dent of and Ti. In the latter regime, time superposition 
holds, i.e. by scaling time with a decay time t'^'^ the Cq 
for the different values of t„ fall onto a master curve at 
a given wavenumber q. This behavior is very similar to 
that found for the KALJ mixture. However, it is differ- 
ent from the behavior predicted by mean-field spin glass 
models and the activated dynamics scaling, as proposed 
by Wahlen and Rieger. Thus, these results suggest that 
the aging dynamics in silica, the prototype of a glass- 
former with a tetrahedral network structure, is very sim- 
ilar to that of simple glass-formers with a closed-packed 
hard-sphere-like structure. We find a difference between 
the KALJ mixture and Si02, however, in the parametric 
plot oiC'^{Cq). For Si02 C'^{Cq) shows a data collapse for 
different sufficiently large and thus Cq — C{q, z{t^, t)) 
whereas this data collapse does not hold as well in the 
case of the KALJ mixture. 

The rest of the paper is organized as follows: In the 
next Sec. we give the details of the BKS potential and 
the simulation. Then, we present the results in Sec. Ill, 
before we summarize in Sec. IV. Appendix A describes 
the implementation of the Nose- Hoover thermostat used 
in our simulation. 



II. MODEL AND DETAILS OF THE 
SIMULATION 

The interactions between the particles are modeled by 
the BKS potential which has been used frequently 
and has proven to be reliable for the study of the dy- 
namics of amorphous silica [25l - l39| The functional form 
of the BKS potential is given by a sum of a Coulomb 
term, an exponential and a van der Waals term. Thus 
the potential between particles i and j, a distance 
apart, is given by 



A: 



C 



(1) 



where e is the charge of an electron and the con- 
stants Aij, B 



and Cy are ^siSi = 0.0 eV, ^siO 



1388.7730eV, ^siSi = 0.0A-\ 
Boo = 2.76000 A-1, Csisi - 
133.5381 eVA-6 and Coo 
The partial charges qi are 



-1.2 and is given by 



18003. 7572eV, Aoo = 
Bsio = 4.87318A-1, 
0.0eVA"6^ Csio = 
175.0000 eVA-6 [4|. 
gsi = 2.4 and go = 
1602.19/(47r8.8542)eVA. 

The Coulombic part of the interaction was computed 
by using the Ewald method [4l|, |4^ with a constant 
aL = 6.3452, where L is the size of the cubic box, and 
by using all g-vectors with |g| < 6 • 2ti/L. [Fi']. We 
ensure that the Ewald term in real space is also differ- 
entiable at the cutoff by smoothing similarly to Eq. (3) 
in 43] with Tc = 8 A and c? = 0.05 A^. To increase com- 
putation speed the non-Coulombic contribution to the 
potential was truncated, smoothed and shifted at a dis- 
tance of 5.5 A. Note that this truncation is not negligible 
since it affects the pressure of the system. In Ref. [i^ 
further slight variations on the potential are described 
in detail ,52] . In order to minimize surface effects pe- 
riodic boundary conditions were used. The masses of 
the Si and O atoms were 28.086 u and 15.9994 u, re- 
spectively. The number of particles was 336, of which 
112 were silica atoms and 224 were oxygen atoms. For 
all simulation runs the size of the cubic box was fixed 
to T = 16.920468 A which corresponds to a density of 
p — 2.323g/cm'^, a value that is very close to the one of 
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real sihca glass, p — 2.2g/cm^ 

We investigated the aging dynamics for systems which 
were quenched from a high temperature T; to a low tem- 
perature Tf . To increase the statistics, for each (Ti, Tf ) 20 
independent simulation runs were performed. To obtain 
20 independent configurations we carried out molecular 
dynamics (MD) simulations using the velocity Verlet al- 
gorithm with a time step of 1.6 fs at 6000 K. The tem- 
perature was kept constant at 6000 K with a stochastic 
heat bath by replacing the velocities of all particles by 
new velocities drawn from the corresponding Boltzmann 
distribution every 150 time steps. Independent configura- 
tions were at least 3.27 ns apart. Each of these configura- 
tions undergoes the following sequence of simulation runs 
(see also Fig. [1]). After fully equilibrating the samples at 
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(Nose Hoover) (Velocily Verlet) 



waiting time tw measurement 

FIG. 1: Schematic sketch of the protocol for the simulation 
runs. 20 independent initial configurations are obtained via one 
long simulation run at 6000 K. For each independent config- 
uration the system is quenched instantaneously and then fully 
equilibrated at temperatures Ti — 3760 K and 5000 K, followed 
by instantaneous quenches to the temperatures Tf = 3250 K, 
3000 K, 2750 K, and 2500 K. At each Tf, time-dependent cor- 
relation functions are determined for different waiting times t^. 
Temperature is kept constant at Tf by coupling the system to a 
Nose-Hoover thermostat. The thermostat is switched off after 
0.33 ns, followed by the continuation of the simulations in the 
microcanonical ensemble for 33 ns. 
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FIG. 2: (Color online) Temperature T = ^^j^'" as a function 
of time t for T = 5000 K, T = 2500 K, 2750 K, 3000 K, 3250 K 
shown in each case for the 11th independent simulation run. 
Skin includes a time average and error bars indicate the corre- 
sponding fluctuations. 



III. RESULTS 



the initial temperatures Tj — 5000 K (for 16.35 ns) and 
Ti = 3760 K (for 32.7 ns), the system was quenched in- 
stantaneously to Tf e {2500 K, 2750 K, 3000 K, 3250 K}. 
To disturb the dynamics minimally, we used a Nose- 
Hoover thermostat [i^ instead of a stochastic heat 
bath to keep the temperature at Tf constant. A velocity 
Verlet algorithm was used to integrate the Nose-Hoover 
equations of motion (see Appendix |X| with a time step 
of 1.02 fs. After 0.33 ns the Nose-Hoover thermostat was 
switched off and the simulation was continued in the NVE 
ensemble for 33 ns using a time step of 1.6 fs. Whereas 
previous simulations used instead the NVT ensemble for 
the whole simulation run, we chose to switch to the NVE 
ensemble to minimize any influence on the dynamics due 
to the chosen heat bath algorithm. For the comparison 
with previous simulations and to check for the lack of 
a temperature drift, we show in Fig. [2] for exemplatory 
simulation runs the temperature T — ^^j^'" as a function 
of time where i?kin is the time averaged kinetic energy 
with fluctuations as indicated with error bars. We find 
that even after switching off the heat bath [s^ there is no 
temperature drift for T = 3250 K and T = 3000 K and for 
T = 2500 K and T = 2750 K there is only a slight temper- 
ature drift which is of the same order as the temperature 
fluctuations and the drift occurs only for t ^ 0.6 ns. For 
all times t ^ 0.6 ns and for all investigated temperatures 
there is no temperature drift and thus the comparison 
with previous simulations valid. 



In all following, we investigate how the structure and 
dynamics of the system depend on the waiting time 
elapsed after the quench from Ti to Tf. We varied the 
waiting time in the range Ons <t^ < 23.98 ns. 

A. Partial Structure Factor 

Figure |3] shows for the temperature quench Tj = 5000 
K to Tf = 2500 K the partial structure factors [ll| 

5.M9,iw) = ^(EEe-(^-(*-)-^^M\ (2) 

\'*=ij=i / 

where r.^ and rj are the positions of particles i and j of 
species a,/? = O, Si. The partial structure factors for all 
other (TijTf) combinations are very similar. Although 
Fig-Elshows S{q,t„) for the largest investigated temper- 
ature quench, there is only a slight t„-dependence for 
very short waiting times < 0.33 ns and almost no t^- 
depcndence for > 0.33 ns. 

B. Generalized Incoherent Intermediate Scattering 
Function 

In this section we focus on the time-dependent gener- 
alized intermediate incoherent scattering function [ll| 
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FIG. 3: (Color online) Partial structure factors Sa()[q,t„) as 
defined in Eq. ([2]) for the temperature quench Ti = 5000 K to 
Tf = 2500 K. Indicated with arrows are the wave vectors q which 
have been used to determine C'q{tvj,tw + t) as defined in Eq. ([B]!. 



which is a measure for the correlations of the positions 
at time t^j and at a later time (tw + t)- We investi- 
gated wave vectors of magnitude q = 1.7, 2.7, 3.4, 4.6, 5.5 
and 6.6 A~^, as indicated with arrows in Fig. [3] We 
show in Fig. m and Figs. [5]- HI results for the first sharp 
diffraction peak at q = 1.7 A~^. Similar results are found 
for all other investigated wave vectors. Figure S] shows 
Gq{tw,tw + t) for the largest investigated temperature 
quench from T, = 5000 K to Tf = 2500 K for waiting 
times iw = - 23.98 ns, as listed in the figure caption of 
Fig. m We find that Cq{tv,,tvj + t) is dependent on 
for all but the last three investigated waiting times. For 
very short times t<5-10^^ns and zero waiting time, 
Cq{tvj = 0, i) is well approximated by Cq of the high tem- 
perature Ti = 5000 K from which the system has been 
quenched (see dashed line in Fig. |4]). Thus, Cg(tw — 0,t) 
for very short times is only dependent on Tj, q and the 
particle type, but independent of Tf. For times of the 
order oi t — 10~'^ns, Cq{t^,t^ + t) is oscillatory due to 
the small system size. For times t > 10"'^ ns, Cq decays 
to zero without forming a plateau for small With 
increasing a plateau is formed, which is independent 
of tw for > 0.33 ns. 

To characterize the plateau height we define F as the 
time average of Cq{tv,,tv! + t) for times 2.55 ps < t < 
6.64 ps. The inset of Fig. [5] shows that for large waiting 
times F{tyj) becomes independent of tw and of Tj. To test 
this independence of Tj further, we show f as a function 
of q for tyj = 16.67 ns in Fig. [Sj We find that the plateau 
height is dependent on the particle type and decreases 
with decreasing g, but F is independent of Tj. 
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FIG. 4: (Color online) Cq{t„,t„ + t) as defined in Eq. ^ for 
the quench Ti = 5000 K to Tf = 2500 K for 0-atoms and q = 
1.7 A^^. Waiting times for the solid lines are from bottom to 
top tw = ns, 1.63 • 10"* ns, 1.63 • 10"^ ns, 1.63 • 10"^ ns, 
0.33 ns, 0.49 ns, 1.17 ns, 1.96 ns, 8.83 ns, 16.67 ns, and 23.98 ns. 
The order of solid lines are for increasing i„ black thin, green 
(gray) thin, black thick, green (gray) thick, black thin, green 
(gray) thin, etc. Error bars are as indicated exemplary. The 
thick dashed line corresponds to Cq{t^,t^ + t) = Fs{q,t) at 
5000 K. 



The plateau in Cq becomes more horizontal with de- 
creasing final temperature Tf, as can be seen by the com- 
parison of Fig.|l(Tf = 2500 K) and Fig.E](Tf = 3000 K). 
Times t > 0.1ns correspond to the a-relaxation, where 
Cq{t„,t^ + t) decays from the plateau to zero. For the 
quench from 5000 K to 2500 K (see Fig. g]) this decay 
depends on for all < 8.83 ns. However, for > 
8.83ns (the largest three iw), Cq{t^,t^+t) is independent 
of not only for intermediate times t (plateau) but for 
all times (including the a- relaxation). Thus, the system 
reaches equilibrium during the simulation run. For the 
quench from 3760 K to 3000 K (see Fig.[6l), Cq{t^,t^ + t) 
becomes independent of for t„ > 1 ns, which means 
that the time at which equilibrium is reached is depen- 
dent on the temperature quench. 

To estimate the time when the system reaches equilib- 
rium for each (Ti,Tf) combination, we next quantify the 
decay of Cq (iw , + • Instead of taking a vertical cut 
in Cq{t^,t^ + t) as we did for T, we now take a hori- 
zontal cut. We define the decay time tjr''^ to be the time 
t = <cq foj. ^jjjgjj Cq{t„,t„ + t^'i) = Cent- We chosc 
Ccut = 0.625/0.41/0.295/0.155/0.085/0.04) for the Si 
particles and C^t = 0.625/0.305/0.195/0.08/0.04/0.014 
for the O particles at the wavenumbers q = 
1.7/2.7/3.4/4.6/5.5/6.61-1, respectively The resulting 
decay times t^'^ as a function of waiting time iw are shown 
in Fig. [7^ for 0-atoms and in Fig.[7}D for Si-atoms. Color 
(black or green/gray) indicates the initial temperature Tj 
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FIG. 5: (Color online) Plateau height, F, as defined in the text, 
as a function of wave vector q (t„ — 16.67 ns) and in the inset as 
a function of tw (? = 1.7 A^^). Error bars are of the same size as 
symbols as indicated exemplary for Ti = 5000 K and Tf = 2500 
K. 
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FIG. 7: (Color online) Decay time t'^'^ for q = 1.7 A"^ and (a) 
0-atoms and (b) Si-atoms. Green (gray) is used for Ti = 3760 K 
and black for Ti = 5000 K. Symbol shape indicates Tf as given 
in the legend. Error bars are given exemplary for (Tj = 3760 K, 
Tt = 2500 K), (Ti = 5000 K, Tt = 2500 K) and (Ti = 5000 K, 
Tt = 3000 K). 



FIG. 6: (Color online) C,(tw, U + t) for the quench Ti = 3760 K 
to Tf = 3000 K. Waiting times and corresponding line styles 
are the same as in Fig. |4l Error bars are of the same order 
as indicated in Fig. O The thick dashed line corresponds to 

C,(tw,tw + t) = fs(g,t) at 3760 K. 



and symbol shape indicates the final temperature Tf. 

We find that <p''(<w) is characterized by three differ- 
ent tw-windows. (I) For waiting times < 0.3 ns, de- 
cay times are significantly lower for Ti = 5000 K (black 
lines and symbols) than for Ti = 3760 K (green/gray lines 
and symbols). The dependence of t^'^ on t^j is strongly 
dependent on all varied parameters, i.e. Tj, Tf, particle 
type, and q. For T; = 5000 K, Tf = 2500 K, 2750 K, and 



q = 1.7 A^^, 2.7 A^^ t^'^{t^) follows roughly a power law 
with an exponent fi « 1.15 with variations of the order 
of 0.07 dependent on Tf, particle type, and q. (II) For 
intermediate waiting times, tp'i(iw) also follows roughly 
a power law with a different exponent than in regime (I) . 
We find for Ti = 5000 K, Tf = 2500 K, 2750 K and q=1.7 
A"\ 2.7 A"i ^ sa 0.35 with variations of the order of 0.08 
depending on Tf, particle type, and q. Best power law 
fits are for Ti = 3760 K, Tf = 2500 K with fi ranging 
from = 0.55/0.57 for q = 1.7 A'^ to ^ = 0.69/0.63 
for q = 6.6 A~^ and for Si/0 atoms. Kob and Barrat 
find for the binary Lennard-Jones system also a power 
law for t^'i(tw)^owever, with n = 0.882 Similar to 
Grigera et al. [43| we find that the transition from small 
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waiting times (I) to intermediate waiting times (II) is 
accompanied by a change of the exponent fj,. (Ill) For 
very long waiting times t^'i(tw) is independent of and 
Tj, i.e. equilibrium is reached. The waiting time ^23 for 
which the transition from regime (II) to regime (III) oc- 
curs is dependent on Tf. ^23 ~ 0.3 ns for Tf — 3250 K, 
^23 « 1 ns for Tf = 3000 K, ^23 « 3 ns for Tf = 2750 K, 
and ^23 w 10 ns for Tf = 2500 K. [Mi] 

Mean-field spin glass models predict d, Q 
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Cq (iw 7 "I" ^) — Cq (^) 



-t) 



(4) 



according to which Cq {t^ + t) can be separated into 
a short-time term C^^ (t) that is independent of and 
an intermediate-time term that scales as h{t^ + t)/h{tvj) 
where h is a monotonously increasing function. It has 
been observed for different systems that the function h{t) 
follows h{t) w t". Thus, the so-called "simple aging" (see 
[igI and references therein) applies and, as a consequence. 
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Miissel and Rieger [481 
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FIG. 8: (Color online) Cgit/t?"^) for the quench from Ti = 

2500 K to Ti = 5000 K for 0-atoms and q = 1.7 Waiting 

times and corresponding lines in the inset are the same as in 

Fig. a 



CqiU,U+t)^C^^{t)+C^'' 



In [(tw + t)/Tfit] 
hr [tw/Tfit] 



(5) 



where the characteristic time scale rgt is a fit parame- 
ter. We find that neither Cq (^7^), nor Cq (t^), nor 

^1 ( '"infl^/rrlf ^ ) (^'^'^ choicc of Tfit) superimposc for 
different t^. Instead we find, similar to the results of Kob 
and Barrat [lB| for a binary Lennard- Jones system, that 
time superposition holds, defined by 



Cq{U,t^+t) = C^'{t)+C^ 



~(AG 



4-Cq/, 



(6) 



In Fig. [SI we show Cq{t/t^'^) for the same set of parame- 
ters as in Fig.m When all waiting times are included (see 
inset of Fig. time superposition does not apply due to 
including too short waiting times. Wahlen and Rieger 
[2(|] have studied Cq(t^,t^ + t) for the same BKS-Si02 
system, however for waiting times smaller than 50 ps. 
Their results are consistent with the inset of Fig. [5] For 
waiting times tw > 0.49ns (see Fig. [5]), however, Eq. © 
is a good approximation. Please note that Cq{t/t'^°') fol- 
lows time superposition for all waiting times > 0.49 ns, 
i.e. not only for the time-range (II), but also for the time- 
range (III). That means for the a-relaxation that the 
shape of the out-of equilibrium curves is the same (within 
error bars) as the shape of the equilibrium curves. We 
find similar results for all other (Ti,Tt) combinations. Si- 
particles, and all other q. 

Next we test whether Cq — C{q,t^,t^ + t) scales as 
Cq = C(g, z{t^, t) where z is a function of and t only, 
i.e. independent of q. Following an approach of Kob and 
Barrat [17[ we show in Figure [HI a parametric plot for 
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FIG. 9: (Color online) C,/ (fw, iw+t) as a function of Cq{t^, t„+ 
t) for q = 1.7 A"^ and q' = 2.7, 3.4, 4.6, 5.5, 6.6 for 0- 
atoms and the temperature quench from 5000 K to 2500 K. 



Cg'(tw,^w + t) as a function of Cq{t„,t„ + t) for q = 
1.7 A-i and q' = 2.7, 3.4, 4.6, 5.5, 6.6 A"! for 0-atoms 
and the temperature quench from 5000 K to 2500 K. For 
sufficiently large iw we find, contrary to the results of 
Kob and Barrat for the Lennard-Jones system, that for 
Si02 the parametric curves superimpose and thus that 
C(g,iw,iw + t) = C{q,z{t^,t) for > 0.49 ns. This 
includes, within error bars, also the equilibrated curves 
for tw }^ 10 ns. We find similar results for all other (Tj, Tf) 
combinations. Si-particles, and all other q. 
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FIG. 10: (Color online) Mean square displacement Ar^(tw, + 
t) as defined in Eq. d?} for the temperature quench from 5000 K 
to 2500 K and for 0-atoms. Waiting times and corresponding 
line styles are the same as in Fig. O 
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FIG. 11: (Color online) tf"^ for 0-atoms. Symbols for the dif- 
ferent (TijTf) combinations are the same as in Fig. [7] Error bars 
are indicated exemplary for (3760 K, 2500 K), (5000 K, 2500 K) 
and (5000 K, 3000 K). 



C. Mean Square Displacement 

In the previous section, we have focused on the analysis 
of Cq(tyj, tvj+t) and identified different time- windows. In 
this section, we consider the mean square displacement 

1 ^ 

i=l 

Figure [TU] shows Ar^(t„,t„ + t) for the temperature 
quench from 5000 K to 2500 K and for 0-atoms. As in 
Fig. m for times t < 5 • 10~^ns and zero waiting time, 
the mean square displacement Ar^(tw = 0,i) is well ap- 
proximated by Ar^ of the high temperature Tj — 5000 K 
from which the system has been quenched (see dashed 
line in Fig. [TU|) and thus independent of Tf . For times 
t « 1Q~^ ns, Ar^{tv,,tvj+t) is oscillatory due to the small 
system size [31}, while for times t > 10"'^ ns and waiting 
times > 0.33 ns, we find that Ar^ forms a plateau 
which is independent of tw As for Cq, we find that the 
plateau is the more horizontal the smaller T{ and the 
plateau height depends on the particle type but is inde- 
pendent of Ti. 

For waiting times > 0.33 ns and times t > 0.1ns, 
the mean square displacement leaves the plateau and in- 
creases further. To characterize the dependence of this 
a-relaxation we define the time i™'^'^ as the time t = t™^'^ 
for which Ar^{t^,t^+tf'^'^) = 1.35 (see Fig.HH). We 
can identify again the three time windows (I) of waiting 
times tv, ^ 0.3 ns with a dependence on Tj, T{ and particle 
type, (II) the aging regime of intermediate waiting times 
where t'^^'^ follows roughly a power law, and (III) for 
very long waiting times when equilibrium is reached. The 
transition from (II) to (III) occurs at approximately the 
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FIG. 12: (Color online) Ar^(f/f™"'^) for the temperature quench 
from 5000 K to 2500 K and for 0-atoms. 



same times t23 as for Cq, i.e. t23 ~ 0.3ns for Tf — 3250K, 
t23 « Ins for Tf = 3000 K, tza « 3ns for Tf = 2750 K 
and t23 ~ 10 ns for Tf = 2500 K. 

Figure [T^] shows the equivalent of Fig. |S] to test time 
superposition. We find for Ar"^ {t / 1^^'^) that time super- 
position is valid for waiting times 0.34 ns < ^ 8.83 ns, 
i.e. for the time window (II) but not for the time window 
(III). 
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IV. SUMMARY 

Using molecular dynamics simulations, we investigated 
for the strong glass former Si02 the aging dynamics be- 
low the critical MCT temperature T^, using the BKS po- 
tential to model the interactions between silicon and oxy- 
gen atoms. After an instantaneous quench from Tj > Tc 
to a temperature T{ < Tc the dynamics towards equilib- 
rium was studied as a function of waiting time tw Note 
that the temperatures Tf were chosen such that equilib- 
rium was reached on the time span of the simulations 
(of the order of 30 ns). The central quantities considered 
in this work are the incoherent intermediate scattering 
function Cq{t^,tw + t) and the mean square displace- 
ment Ar^(tw,tw + t). These functions depend on the 
time origin at as long as is smaller than the typical 
relaxation time, Teq, that is required to equilibrate the 
system. 

We find that the decay of Cq{t^, + 1) (and similarly 
the rise of Ar'^{t^,t^ + t)) exhibit qualitative changes 
from short to long waiting times. At short waiting times, 
relaxation processes are dominant that correspond to the 
early (3 relaxation regime at the target temperature Tf. 
In this regime, no well-defined plateau is found in 
Cq{t^,t^ + t) (see Fig. 2]). Instead, this function first 
decreases rapidly, followed by a strongly stretched expo- 
nential decay to zero. At long waiting times, the /3 re- 
laxation seems to be very similar to that at equilibrium. 
The Debye- Waller factor (i.e. the height of the plateau 
in Cq(iw,tw + t)) has reached its equilibrium value, al- 
though the decay of Cg(tw, ^w+t) from the plateau to zero 
is faster than that at equilibrium. However, the shape of 
curves describing the long-time decay of Cq{t^, t^, + 1) is 
the same as that at equilibrium. Thus, Cq follows a sim- 
ple time superposition for long waiting times, different 
e.g. from the "activated dynamics scaling" proposed by 
Wahlen and Rieger for BKS silica. 

Our results show that the aging dynamics of BKS sil- 
ica is very similar to that of the KALJ mixture. For 
both silica and the KALJ mixture three t„-regimes can 
be identified and Cq follows time superposition for suffi- 
ciently large t^. The only difference between these two 
systems is that Cq scales as C{q, z{t„,t) for Si02 but 
less well for the KALJ mixture. So slightly below its 
critical temperature Tc of MCT, the strong glass- former 
silica does not seem to be very different from typical frag- 
ile systems, although one has already reached the low 
temperature Arrhenius regime (note that activation en- 
ergies for the self-diffusion, viscosity etc. are of the order 
of 5 eV, similar to the corresponding activation energies 
close to the glass transition temperature Tg w 1450 K, as 
measured in various experiments). However, the dynam- 
ics could be very different at very low temperatures (close 
to Tg) where the long-time aging regime is not accessi- 
ble by computer simulations. Thus, more experimental 
work on the aging dynamics of silica around Tg would 
be very desirable. We also leave for future work to test 
whether also for other systems three tw-regimes are iden- 



tified and whether the equilibrium curve is included in 
the time superposition of Cq at intermediate and large 
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Appendix A: Velocity Verlet Nose Hoover 

The Nose Hoover equations of motion are for particles 
i = 1, . . . , N at position ri with momentum pi 

Pi 



(Al) 
(A2) 



and thus 



= Ci"« 

mi 



(A3) 



(A4) 



d^ In J 



= ^ = ^ (^5]m,r2-XA;T^ (A5) 

where X = 3N . We integrated with the generalized 
velocity Verlet form of Fox and Andersen [i^ 



2 \ m 
In s{t + At) = In s{t) + At ({t) 



-muit) 



{Atf ' " 
2Q 



^m,x\{t) ~XkT 



1=1 



CPP™^(t-hAi) = Q(t) 



(A6) 



(A7) 



(A8) 



v^(t + At) = i,(t) (A9) 
^ At /E\(t)_+I\(t^HAt) 
2 V 

-[C(t)+C"(t + At)]r,(t) 

1 _ ^^approx^^ _^ 



At / 



N 



C(t + At) - C(i) + ^(E™»^?W (AlO) 



^ TO.r^ (t + At) -~2XkT) 
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1\ 2 

To ensure that T = ^ ^ + J7({rJ) + + XkTlns 

i=l 

is conserved (see [1^) we chose Q — 50000 u. 
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